clear
rng default
addpath /Users/tew647/mosek/9.3/toolbox/r2015aom
%% Load useful matrices
load PYX_cond 
PYX_cond_split=[PYX_cond(1:3); PYX_cond(4:end)];
load U_CS

%% Construct grid of parameters to check for each realisation of X
gran=0.01;
u0grid_temp=0; %[Ychosen=0, X=1] [Ychosen=0, X=2], location normalisation 
u1grid_temp_all=[U_11_CS U_21_CS]; %[Ychosen=1, X=1] [Ychosen=1, X=2], scale normalisation
u2grid_temp=-20:gran:20;  %[Ychosen=2, X=1] [Ychosen=2, X=2]
% try different values of gran and different sizes of grid hypercube as discussed in Appendix C of the paper (need powerful cluster to run the code for higher values of gran and larger sizes of hypercube)

%% Construct identified set for each realisation of X
ntypes_m=2;
IdSetM=cell(ntypes_m,1);

for x=1:ntypes_m %we can proceed separately across types as long as we do not impose independence from X
    u1grid_temp=u1grid_temp_all(x);
    [ca, cb, cc] = ndgrid(u0grid_temp, u1grid_temp, u2grid_temp);
    u0grid=ca(:);
    u1grid=cb(:);
    u2grid=cc(:);
    sg=size(u0grid,1);
    % Components invariant across parameter values
    [sgb, ...
    equalorder, ...
    Ugridreduced, sgreduced,...
    id_Ugridreduced,...
    s_id_gridreduced, sU_id_gridreduced,...
    n_U_sets, Tcoord, indices_pairs,...
    numeqconstraints,numzeros,coordrowseq,d3, e5, f5, g3, h5, i5, l3,m5, n5, fillAeq]=useful_anyparam3(u0grid, u1grid, u2grid, sg);  
    % Loop over parameter values, order of columns: [u0x u1x u2x]
    IdSetM{x}=OneSide3(PYX_cond_split(x,:), u0grid, u1grid, u2grid,...
                        sgb, ...
                        equalorder, ...
                        Ugridreduced, sgreduced,...
                        id_Ugridreduced,...
                        s_id_gridreduced, sU_id_gridreduced,...
                        n_U_sets, Tcoord, indices_pairs,...
                        numeqconstraints,numzeros,coordrowseq,d3, e5, f5, g3, h5, i5, l3,m5, n5, fillAeq);  
end

[cx1, cx2] = ndgrid(1:size(IdSetM{1},1), 1:size(IdSetM{2},1));
cx1=cx1(:);
cx2=cx2(:);
IdSetM_final=[IdSetM{1}(cx1,:) IdSetM{2}(cx2,:)]; 
save('IdSetM_final.mat', 'IdSetM_final')


